function preset()
{
    #=================================
    fig_name=3d_line_drawing
    fig_fmt=png
    #=================================
    bathy=../../gmt/Geersen_2018_Tectonics.nc
    bathy_merge=../../gmt/P1_3D_merge.grd
    aftershock=../interface/flo_aftershock/after_latlongdep.txt
    main=../interface/Schurr2014/Schurr_2014_Nature_main.grd
    after=../interface/Schurr2014/Schurr_2014_Nature_after.grd
    cdp_num=cmp_num_25.txt
    # 透明度
    alpha_profile=50
    #图片大小
    fig_x=8
    fig_y=8
    fig_z=2
    #数据范围
    zmin=-40000
    zmax=0
    lon_min=-71.9
    lon_max=-70.2
    lat_min=-21.3
    lat_max=-18.8
    # 视角
    ##方位角
    azimuth=150
    ##仰角
    elevation=15
    view=$azimuth/$elevation
    #range
    range_lon_lat=$lon_min/$lon_max/$lat_min/$lat_max
    zrange=$zmin/$zmax
    range=$range_lon_lat/$zrange
    
}

preset

gmt begin $fig_name $fig_fmt E1000 TG
    #==========zmin作图==================
    #gmt grdcut $bathy_merge -R$range_lon_lat -V -Ggrdcut.grd
    #gmt grdgradient grdcut.grd -Nt1 -A60 -Gintens_merge.grd -V
    gmt basemap -R$range -JM${fig_x}c -JZ${fig_z}c -Ba -BwSEn -pz$view/$zmin \
    --FONT_LABEL=8p,Helvetica-Bold --MAP_ANNOT_OFFSET_PRIMARY=0.05c --FONT_ANNOT_PRIMARY=8p,Helvetica-Bold --MAP_FRAME_PEN=0.5p,black #3D_framework # -Bza20000f2000+l"Depth(m)" -BwSEnZ+b
    depth_km=`echo "scale=2; $zmin / 1000 " | bc`
    #gmt makecpt -T-10000/1000/5000+ -Cgray > bathy_gray.cpt
    gmt grdimage $bathy_merge -R$range -JM${fig_x}c -JZ${fig_z}c -pz$view/$zmin -Cbathy_gray.cpt -Iintens_merge.grd -Bwsen  #bathymetry
    gmt basemap -R${range_lon_lat}/${depth_km}/0 -JM${fig_x}c -JZ${fig_z}c -Bza10f1+l"Depth(km)" -BwSEnZ+b -pz$view \
    --FONT_LABEL=8p,Helvetica-Bold --MAP_ANNOT_OFFSET_PRIMARY=0.05c --FONT_ANNOT_PRIMARY=8p,Helvetica-Bold --MAP_FRAME_PEN=0.5p,black #-BZ

    for v in 04 25 15 #dipping line and strike line
    do
        gmt plot ../../gmt/P1_basemap_shade/MGL1610MC${v}_new.input -R$range -JM${fig_x}c -JZ${fig_z}c -pz$view/$zmin -W0.04c,black 
    done
    for v in 04 25 15  #cmp num and cmp marks
    do 
        awk '{if($1>=3000 && $1<=25000) print $0}' ../../gmt/P1_basemap_shade/MGL1610MC${v}.latlong > tmp.txt
        awk 'NR%2000==1{print $2,$3,$1}' tmp.txt  > tmp1.txt
        gmt text tmp1.txt -R$range -JM${fig_x}c -JZ${fig_z}c -pz$view/$zmin -F+a-90+f5p,Helvetica-Narrow-Bold -D0c/-0.3c
        gmt plot tmp1.txt -R$range -JM${fig_x}c -JZ${fig_z}c -pz$view/$zmin -W0.03c,red -Sc-0.05c  
    done
    gmt coast -W0.25p,black -Ggrey -JM${fig_x}c -JZ${fig_z}c -pz$view/$zmin

    #==========空间作图==================
    for v in 04 25 15  #seafloor
    do 
        gmt plot3d ../slice_seismic/seafloor_val_${v}.txt -Sp0.1p -R$range -JM${fig_x}c -JZ${fig_z}c -pz$view #plot seafloor 
        gmt plot3d ../slice_seismic/seafloor_vertical_start_${v}.txt -W0.1p,black -R$range -JM${fig_x}c -JZ${fig_z}c -pz$view # draw vertical lines at beginning
        gmt plot3d ../slice_seismic/seafloor_vertical_end_${v}.txt -W0.1p,black -R$range -JM${fig_x}c -JZ${fig_z}c -pz$view #draw vertical lines at end
    done
    for v in 25 04 #horizon
    do
        awk '$3~/[[:digit:]+]/{print $1,$2,$4*-1,$3}' ../line_drawing/${v}_horizon.dat > tmp1.txt #horizon
        ZONE=19
        MERIDIAN=`echo $ZONE | awk '{print int((('$ZONE'*6)-180)-3)}'`
        UTM_REGION=`echo $MERIDIAN| awk '{printf("%i/%i/%i/%i\n", $1-3,$1+3,-90,0)}' `
        awk '{if(NR>=1) printf("%.12lg %.12lg %.5f %d \n",$1,$2,$3,$4)}' tmp1.txt | \
        gmt mapproject -V -Ju$ZONE/1:1 -R$UTM_REGION -C0/0 -F -I | \
        awk '{printf("%.12lg %.12lg %.5f %d \n",$1,$2,$3,$4)}' > tmp0.txt #转换 utm2latlong
        k=$(awk 'END{print $4}' tmp0.txt) #the max num of the layers
        i=$(awk 'NR==1{print $4}' tmp0.txt)
        while [ ${i} -le ${k} ]; do
            awk '$4 == '${i}'{print $1,$2,$3}' tmp0.txt > tmp.txt
            gmt plot3d tmp.txt -W0.1p,magenta -R$range -JM${fig_x}c -JZ${fig_z}c -pz$view
            i=$(($i + 1))
        done 
    done
    for v in 15 #horizon,因为 15 有两段
    do
        awk '$3~/[[:digit:]+]/{print $1,$2,$4*-1,$3}' ../line_drawing/${v}_horizon1.dat > tmp1.txt #basement
        ZONE=19
        MERIDIAN=`echo $ZONE | awk '{print int((('$ZONE'*6)-180)-3)}'`
        UTM_REGION=`echo $MERIDIAN| awk '{printf("%i/%i/%i/%i\n", $1-3,$1+3,-90,0)}' `
        awk '{if(NR>=1) printf("%.12lg %.12lg %.5f %d \n",$1,$2,$3,$4)}' tmp1.txt | \
        gmt mapproject -V -Ju$ZONE/1:1 -R$UTM_REGION -C0/0 -F -I | \
        awk '{printf("%.12lg %.12lg %.5f %d \n",$1,$2,$3,$4)}' > tmp0.txt #转换 utm2latlong
        k=$(awk 'END{print $4}' tmp0.txt) #the max num of the layers
        i=$(awk 'NR==1{print $4}' tmp0.txt)
        while [ ${i} -le ${k} ]; do
            awk '$4 == '${i}'{print $1,$2,$3}' tmp0.txt > tmp.txt
            gmt plot3d tmp.txt -W0.1p,blue -R$range -JM${fig_x}c -JZ${fig_z}c -pz$view
            i=$(($i + 1))
        done 
        awk '$3~/[[:digit:]+]/{print $1,$2,$4*-1,$3}' ../line_drawing/${v}_horizon2.dat > tmp1.txt #basement
        ZONE=19
        MERIDIAN=`echo $ZONE | awk '{print int((('$ZONE'*6)-180)-3)}'`
        UTM_REGION=`echo $MERIDIAN| awk '{printf("%i/%i/%i/%i\n", $1-3,$1+3,-90,0)}' `
        awk '{if(NR>=1) printf("%.12lg %.12lg %.5f %d \n",$1,$2,$3,$4)}' tmp1.txt | \
        gmt mapproject -V -Ju$ZONE/1:1 -R$UTM_REGION -C0/0 -F -I | \
        awk '{printf("%.12lg %.12lg %.5f %d \n",$1,$2,$3,$4)}' > tmp0.txt #转换 utm2latlong
        k=$(awk 'END{print $4}' tmp0.txt) #the max num of the layers
        i=$(awk 'NR==1{print $4}' tmp0.txt)
        while [ ${i} -le ${k} ]; do
            awk '$4 == '${i}'{print $1,$2,$3}' tmp0.txt > tmp.txt
            gmt plot3d tmp.txt -W0.1p,blue -R$range -JM${fig_x}c -JZ${fig_z}c -pz$view
            i=$(($i + 1))
        done 
    done
    for v in 25 04 #basement
    do
        awk '$3~/[[:digit:]+]/{print $1,$2,$4*-1,$3}' ../line_drawing/${v}_basement.dat > tmp1.txt 
        ZONE=19
        MERIDIAN=`echo $ZONE | awk '{print int((('$ZONE'*6)-180)-3)}'`
        UTM_REGION=`echo $MERIDIAN| awk '{printf("%i/%i/%i/%i\n", $1-3,$1+3,-90,0)}' `
        awk '{if(NR>=1) printf("%.12lg %.12lg %.5f %d \n",$1,$2,$3,$4)}' tmp1.txt | \
        gmt mapproject -V -Ju$ZONE/1:1 -R$UTM_REGION -C0/0 -F -I | \
        awk '{printf("%.12lg %.12lg %.5f %d \n",$1,$2,$3,$4)}' > tmp0.txt #转换 utm2latlong
        k=$(awk 'END{print $4}' tmp0.txt) #the max num of the layers
        i=$(awk 'NR==1{print $4}' tmp0.txt)
        while [ ${i} -le ${k} ]; do
            awk '$4 == '${i}'{print $1,$2,$3}' tmp0.txt > tmp.txt
            gmt plot3d tmp.txt -W0.1p,black -R$range -JM${fig_x}c -JZ${fig_z}c -pz$view
            i=$(($i + 1))
        done 
    done
    for v in 25 04 #fault
    do
        awk '$3~/[[:digit:]+]/{print $1,$2,$4*-1,$3}' ../line_drawing/${v}_f.dat > tmp1.txt 
        ZONE=19
        MERIDIAN=`echo $ZONE | awk '{print int((('$ZONE'*6)-180)-3)}'`
        UTM_REGION=`echo $MERIDIAN| awk '{printf("%i/%i/%i/%i\n", $1-3,$1+3,-90,0)}' `
        awk '{if(NR>=1) printf("%.12lg %.12lg %.5f %d \n",$1,$2,$3,$4)}' tmp1.txt | \
        gmt mapproject -V -Ju$ZONE/1:1 -R$UTM_REGION -C0/0 -F -I | \
        awk '{printf("%.12lg %.12lg %.5f %d \n",$1,$2,$3,$4)}' > tmp0.txt #转换 utm2latlong
        k=$(awk 'END{print $4}' tmp0.txt) #the max num of the layers
        i=$(awk 'NR==1{print $4}' tmp0.txt)
        while [ ${i} -le ${k} ]; do
            awk '$4 == '${i}'{print $1,$2,$3}' tmp0.txt > tmp.txt
            gmt plot3d tmp.txt -W0.1p,red -R$range -JM${fig_x}c -JZ${fig_z}c -pz$view
            i=$(($i + 1))
        done 
    done
    for v in 25 04 #rotat
    do
        awk '$3~/[[:digit:]+]/{print $1,$2,$4*-1,$3}' ../line_drawing/${v}_rotat.dat > tmp1.txt 
        ZONE=19
        MERIDIAN=`echo $ZONE | awk '{print int((('$ZONE'*6)-180)-3)}'`
        UTM_REGION=`echo $MERIDIAN| awk '{printf("%i/%i/%i/%i\n", $1-3,$1+3,-90,0)}' `
        awk '{if(NR>=1) printf("%.12lg %.12lg %.5f %d \n",$1,$2,$3,$4)}' tmp1.txt | \
        gmt mapproject -V -Ju$ZONE/1:1 -R$UTM_REGION -C0/0 -F -I | \
        awk '{printf("%.12lg %.12lg %.5f %d \n",$1,$2,$3,$4)}' > tmp0.txt #转换 utm2latlong
        k=$(awk 'END{print $4}' tmp0.txt) #the max num of the layers
        i=$(awk 'NR==1{print $4}' tmp0.txt)
        while [ ${i} -le ${k} ]; do
            awk '$4 == '${i}'{print $1,$2,$3}' tmp0.txt > tmp.txt
            gmt plot3d tmp.txt -W0.1p,gray -R$range -JM${fig_x}c -JZ${fig_z}c -pz$view
            i=$(($i + 1))
        done 
    done

gmt end show
#rm tmp*